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("^ , We develop renormalization group (RG) methods for solving partial and stochastic differential 

^S| . equations on coarse meshes. RG transformations are used to calculate the precise effect of small scale 

dynamics on the dynamics at the mesh size. The fixed point of these transformations yields a perfect 

operator — an exact representation of physical observables on the mesh scale with minimal lattice 

ry^ • artifacts. We apply the formalism to simple nonlinear models of critical dynamics, and show how 

the method leads to an improvement in the computational performance of Monte Carlo methods. 
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I . I. INTRODUCTION 

"^ ' The purpose of this paper is to introduce numerical methods that avoid unnecessary discretization — or, over- 

1 t ] discretization — purely for the purpose of obtaining adequate accuracy. An important and classical example of this 
Ci • is large eddy simulation in the modeling of turbulent flows. Many large scale flows of engineering, geophysical or 
H I atmospheric interest contain many length scales down to the dissipation scale, yet it is large scale drag that one wants 

I to compute. In such a situation, it is wasteful and undesirable to expend computer time on details that are of no 
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intrinsic interest. 

The approach outlined in this paper builds upon our previous work y to use RG methods to integrate out the 
fj ' dynamics one wishes to ignore, so that numerical methods can instead focus on the appropriate scale of interest. This 
—''^ is not trivial because of scale interference: the nonlinear amplification of the effect of small scale dynamics, which 
contaminates and eventually pollutes the large scale dynamics. There are several distinct facets to this problem. 
First is the representation of the small scale dynamics as a stochastic field that acts on the coarse-grained degrees of 
Qs 1 freedom. As discussed in our earlier paper, this inevitably leads to non-locality. We will see here that it is possible not 
T^ ■ only to coarse-grain individual operators, as in rcf. ||l[ but also to coarse-grain at the level of the governing differential 
^^ , equation. This leads to a theory that is non-local in space and time. This applies to systems with a finite number of 
On ■ degrees of freedom, as well as spatially extended systems, which are the main focus of our work here. 
^^ ' Second, the representation of the theory on the lattice can be improved by systematically integrating out the 

small scales, leading to an effective theory that has no (or few) residual discretization artifacts. This is referred to 
as a "perfect theory" in the literature. We demonstrate how this arises and exhibit this feature by calculating the 
"jrt ■ dispersion relation of the effective theory in the perfect representation. 

r^ ' Our work is related to that of Chorin and co-workers y^jOJI who use optimal prediction methods to treat the lack 

i' . of resolution of small scales. The main differences are that they assume that the small scales are initially in thermal 
'■^ ' equilibrium, and also that they do not attempt to remove lattice artifacts. There has also been an attempt to use 
^ [ similar methods in the study of isotropic turbulence H . 

O • Our work grew out of attempts to improve lattice gauge theory, pioneered by the paper of Hasenfratz and Nieder- 

meyer. For a review of this body of work, the reader is referred to the review article by Hasenfratz S. In addition to 
the work in ref. [|l|, there have been two attempts ||,|l^ to solve differential equations using perfect operators. As we 
will see below, it is not enough to perfectly coarse-grain the individual operators appearing in a partial differential 
equation: once there is a non-inflnitesimal time step, coarse-graining introduces memory effects, so that the entire 
C^ ' differential equation must be represented as coarse-grained in space-time. In addition, it should always be remem- 
bered that there is no unique perfect operator for a given differential operator. A specification must be made of the 
microscopic probability distribution for the small-scale degrees of freedom. These papers implicitly impose a Gaussian 
free field theory distribution on the small-scale degrees of freedom. The methods given in the present article are more 
general, and make no such assumption, explicit or implicit. 

Let us now introduce the problem of removing lattice artifacts. Suppose the dynamics of a spatially extended 
system is described by a partial differential equation (PDF), which yields the solution u(x,t). A standard procedure 
is to sample u{x, t) at points Xi, tj, which are equidistant with spacings Ax and At, and find a discretized form of the 
PDE that is devised to approximate the values Uij = u{xi,tj). The requirement is that in the continuum limit, the 
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sequence Uij converges to u{x,t). The conventional way of discretizing the PDE is to approximate differentiations 
with finite differences. 

The disadvantage of this uniform sampling (US) approach is that one is forced to reproduce as faithfuUy as possible 
all the detail and fine structure of the solution, even on a scale that may be of no interest or worse, beyond the regime 
of applicability of the differential equation itself. This has two consequences: 

• a small grid size Aa; must be used, which implies many grid points must be calculated and stored; 

• for dynamic problems, a small time step At is implied by the small Ax, either for reasons of accuracy or stability 
of the numerical method. 

As a result, there is a huge computational cost associated with this conventional numerical scheme, which makes 
the study of problems such as critical dynamics and pattern formation very difficult to carry out. There is a need for 
improved, physically motivated methods for numerical experiments. 

The purpose of a numerical simulation is to study the macroscopic properties of a physical system. Different 
microscopic dynamics may be related, via coarse graining (CG), to the same macroscopic dynamics that defines a 
universality class. Often CG means the local averaging of a continuous variable, 

/•A/2 

U{X)^ / dxu{X + x), (1) 

where u{x) is the continuum variable, U(X) is its coarse-grained counterpart and A is the coarse graining length scale. 
Instead of focusing on the small scale degrees of freedom, we should determine and use the coarse grained description 
of the system appropriate at the macroscopic scale. 

One of these physics-motivated numerical methods is the cell dynamical scheme |Il| , in which a discrete description 
of the system dynamics is obtained directly from considerations of the underlying symmetry and conservation laws. 
It has been successfully used to tackle problems such as asymptotic scaling behavior in spinodal decomposition Il2] 
and the approach to equilibrium in systems with continuous symmetries, such as XY magnets [13| and liquid crystals 
WM. There have also been attempts at using the RG in dynamic Monte Carlo simulations |p5|, f6|. 

To investigate what is required to obtain a coarse-grained dynamic description, suppose that we denote the coarse- 
graining operator at scale A by the symbol Ca, which transforms u{x,t) to U{X,t). Then conceptually we need to 
find the operator L\ which connects U{X, 0) and U{X, t) given the microscopic time evolution operator L connecting 
m(x,0) with u{x,t), as shown schematically in the commutativity diagram below: 

u{x,0) — > u(x,t) 



Ca 



Ca 



Ua{X,0) ^ UA{X,t) 

Notice that there is not a unique choice of Ca. The usual choice is local averaging. In principle, other operators 
can be used, such as the majority rule scheme used in the coarse graining of Ising spins in thermal equilibrium. Once 
a coarse-graining operator Ca has been defined, there should be a unique prescription to obtain La ^^- In this 
paper, coarse graining is understood to mean local averaging. Later, stochastic coarse graining will be introduced as a 
variant of the simple local averaging. In the development of the theory of perfect operators a parameter, denoted here 
by kq (see section IIIB), naturally arises, which characterizes the nature of the coarse-graining procedure. The form 
(|l|) is only appropriate if kq is infinite; if kq < oo, additional noise terms are generated which reflect the reduction 
in the number of degrees of freedom in the system. As already stressed in our earlier paper IpJ, it is inconsistent to 
work with a perfect operator with kq < oo and to use the kq — oo form Im) as some authors P,EO[ have done. We also 
see no reason why coarse-grained equations should be derived by varying a coarse-grained action in the absense of a 
small parameter, which is the starting point of these authors. Instead we begin with a dynamics which is intrinsically 
stochastic and study the effect of CG on this system. The well-known path-integral formulation of such equations 
may then be used to carry out the CG: there is no need to invoke a variational principle. 

We need to consider the appropriate coarse-graining scale. Two situations are possible here. In the first, we suppose 
that the solution we wish to obtain has a natural scale A below which there is no significant structure. In that case, 
our goal is to avoid having to over-discretize the problem merely in order to attain the accuracy of the continuum limit. 
Thus, we would like to be able to use as large a value for the grid spacing Ax as possible without sacrificing accuracy. 
In the second situation, there is no such obvious scale, or at least, it is not known a priori, but the computational 



demands are so large that it is simply not feasible to work with a grid spacing Ax smaller than some size A. In this 
case, we would like to minimize in some sense the artifacts that must inevitably arise. 

The first situation is more straightforward because the only issue is speed of convergence to the continuum limit: 
there is no explicit discarding of important dynamical information. In the second situation, one is making an un- 
controlled and potentially severe truncation of the correct dynamics. One has to ask: can one model the neglected 
unresolved scales as effective renormalizations of the coefficients in the original PDE? Are the neglected degrees of 
freedom usefully thought of as noise for the retained large-scale degrees of freedom? And how can any available sta- 
tistical information on the small-scale degrees of freedom be used to improve the numerical solution for the large-scale 
degrees of freedom? 

The plan of the paper is as follows. In section II we set up the coarse-graining algebra, which forms the basis of our 
approach, using the path-integral formulation of stochastic dynamics as our starting point. This formalism is then 
used in section III to obtain the perfect operator for dynamics governed by linear operators. Section IV describes the 
results of numerical simulations using the perfect operator with Langevin dynamics and section V using the Monte 
Carlo approach. A range of issues is discussed, from applications of the method to the diffusion equation and nonlinear 
model A dynamics to the question of the truncation of perfect operators required when carrying out simulations. Our 
conclusions are presented in section VI and the structure of the coarse-graining algebra is discussed in an appendix. 

II. COARSE GRAINING IN THE PATH INTEGRAL FORMULATION OF LANGEVIN DYNAMICS 

In this section, we derive the path integral formulation of the Langevin dynamics and present the general framework 
under which the perfect linear operator is derived. The analysis is applicable to both PDEs and stochastic differential 
equations. For simplicity, we study a system whose dynamics is described by a stochastic differential equation (SDE) 
with the following form, 

^^|i^ = -/(x,t;{0}) + 77(x,i), (2) 

where </> is a field, / is the forcing term (it can depend on (p and/or its spatial derivatives), and 77 is a white noise. 

It is convenient to regularize the problem on a (fine) N x N' lattice with grid size Ax and At in the space and time 
directions, respectively. In the lattice picture, all variables in the original PDE are vectors of functions of discrete 
space X — iAx and time t — jAt where i £ [0,N — 1], j G [0, A^' — 1]. We define g(iAx,jAt) = g{i,j) and denote the 
space-time volume element AxAt by AV. The noise satisfies {r]{i,j)) — and (rjiiTJ) ri{i' ,j')) — -^ 6ij Siiji where fl 
is the noise strength and Si^ is the Kronecker symbol. Given the system is in the state 0o a-t time tg, the probability 
that the system will be in state ^1 at time ii is given by p8[ , 

P{(buh\ct>o,to) = JDcpDrj exp{-^ ^ [772(1,^) - ^d^f]} Sirj - dt(b - fm , (3) 

where the integration is over all configurations beginning at ^o and ending at ^i. We can use this path-integral 
formula to determine the dynamics followed by the coarse grained or uniform sampled variable. 

By a discretization scheme, we will mean a process made up of a series of magnifying operations which lead from a 
microscopic description of a system to a macroscopic description on a lattice. These magnification operations are, by 
default, magnification of a length scale by a factor of 2. Coarse graining and uniform sampling are both special cases 
of a discretization scheme. 

Suppose a system is specified by the values of a function /, such as a field configuration, on a fine lattice with 2N 
grid points x — {xi,X2, ■ ■ ■ ,X2n) separated by grid size Ax. One step (level) of coarse graining is defined as local 
averaging of the function's values at every two neighboring sites. 

fn = ^(/2n-l + f2n) , fn = -^{hn ~ /2n-l)- (4) 

Vector / is the coarse grained version of /, while / stores the detailed information that is lost after coarse graining. 
After one level of CG, the system is described by a new function / on a coarser lattice with N grid points separated by 
twice the original grid size of Ax*^ = 2Aa:, where the superscript M indicates "magnified value" . We define 2A^ x N 

projection matrices i?, R such that 



I f = R f , f = R /• 

These matrices act as projection and inverse projection operators between the original functional space and the coarse 
grained functional space. They facilitate an easier mathematical formulation. Many of the properties of the matrices 
can be found in the appendix. If we are interested in an operator O on the original grid, then it is possible to define 
four corresponding operators on the coarse-grained grid, which we denote by Oa, Ob, Oc and Od- For instance, 

Od^ R OR. 

The analogous definitions of Oa , Ob and Oc are given in the appendix. 

A similar algebraic scheme can be defined for the uniform sampling transformation, where the projection operator 
samples every other point and discards the rest: 

fn = /2ri-l, fn = f2n- (6) 

Rrn,n = Srn,2n-1, Rm,n ^ 5m,2n, m (^ [l,2N], ne[l,7V]. (7) 

Using the notations listed above, we can write down the magnification procedure in space for the 1 + 1 dimensional 
version of (g), coarse graining in space only. The integrations over the (p and r] variables are decomposed into 
integrations over </>, (j), fj and fj variables and the rj integration carried out using the delta-function. The remaining 
delta-function is replaced using the identity S{x) = a S{a x) = a J dqe^'' '^^/2tt. This leads to a path integral, neglecting 
any constant factors, of the form 



/AV '2^ 
D^Df^Dq cxp{ ^^ E [- ^' - *<? (^ - ^*^)]} ^ 

y i?^exp{-^^ ^[-(a,<^ + /f + zg/-^^(V + V)]}' (8) 



where the constant c is 1 or 2 for CG and US respectively due to their different projection matrix properties and 
where AV^ = 2AxAt — 2AV^ is the magnified volume element. The important point is that, in general, both / and 
/ are functions of cf) and 4>. 

What we would like to do, is integrate over the </> degrees of freedom, carry out the q integration and end up with a 
form similar to the one we started with, but with new, renormalized, parameters. More specifically, we would like the 
integration over (p to give a result of the form exp{— -^^(igF — ^^aj dj,F)}. Then we could readily integrate over 
q and compare the result with the path-integral form to read off the evolution equation for the new coarse grained 
variable as dt(f> = —F{(j)) + fj. However, we would not expect to be able to do this in general, and as usual in all 
applications of the RG, an approximation scheme has to be developed alongside this formalism in order to make any 
progress. There is, however, one case in which the integrations can be carried out, and that is the linear case. We 
therefore study this first, before returning to the nonlinear case later. 

III. PERFECT OPERATOR FOR DYNAMICS 

In this section, we will determine perfect operators of dynamics governed by linear operators. We will find the fixed 
point fiow of operators for the diffusion equation under CG and US transformations. In addition, the perfect operator 
in discrete space and time is obtained for the diffusion equation and its properties discussed. 

A. Iterative Relations and Fixed Points in the linear case 

We begin by performing the magnifying transformation on the SDE ^) where / is a linear function of (p, that is, 

dt4>^-L(t> + r,, (9) 



where 77 is a white noise. Here _£ is a general hnear operator and contains spatial, but not temporal, derivatives. 
It is assumed to possess inversion symmetry and translational invariance. For the diffusion equation, L is the finite 
difference Laplacian operator with a minus sign. The conventional choice is the central difference operator L™ „ — 

(2 5m,n — 5m^n+l — '5m,n-l)/Ax . 

To obtain the dynamics of the coarse grained variable, we have to integrate out the small length scale degrees of 
freedom in equation (f^. In the linear case, the Jacobian term is constant and so does not enter into the analysis. 
Applying the projection matrices to equation (S), inserting the result into the path integral in equation (Rh and 
integrating out the 4> and q degrees of freedom yields 

/AT/M |- - 

D^Z?^exp{-^^^[f + (^-r;^^)Q-i(^-r?^)J}, (10) 

where 77*^ = dtcf) + {La - LcM~^Ld)(I> and Q = LcM^^iM'^y^L'^. Here the operator M is given by Idt + Lb- 
Defining a new noise source ff = fj — [I + Q]~^ri'^^ and carrying out the integration over fj' yields 

/AyM 
D^Dij'' exp{-^ ^ 77*^(7 + Q)-W} S{v^' - dt^ - {La ~ LcM-^Ld)4>) ■ (H) 

Comparing with the form (||) , it follows that the dynamic equation satisfied by 4> is 

where L'-^'~^ = La — LcM^^Ld- The new noise source 77*^ is no longer a white noise: it has a spatial correlation as 
well as a time correlation: 

(77^0-0 and(77*^(r,i)77^^(/,0) = ^(/ + g)(r-r',t-0- (12) 

Given that the noise source is no longer Markovian after the first step of coarse graining, we need to start with a 
more general noise source in order to iterate the coarse graining procedure. Define a general Gaussian noise source 
with the following properties, 

(77) =0 and (7;(r,i)^(r',t')) - -^p'\r - r' ,t - t') . (13) 

Repeating the above analysis, we find that the coarse grained dynamic equation remains the same, however the coarse 
grained correlation matrix is modified and is given by 

(pCG)-i ^ LcM-'p-B\M-'fLl + t{pA^ PcPS'PDr' f ^ , (14) 

where F = / + LcM~^pg pn. The presence of time derivatives in p makes the noise non-Markovian. In general, we 
should be careful about the boundary term in this case [191. In particular, we need to specify corresponding initial 
conditions for each time derivative generated through the iterative relation. 

The first term in L — La — LcM^^Ld is not what we would naively choose as the Laplacian operator with 
a coarse grained grid size Aa;^^. Instead, the second term, which comes from accounting for the influence of the 
integrated out small length scale degrees of freedom, gives an important contribution to the coarse grained operator 
and cannot be treated as a perturbation. 

It is more convenient to examine the coarse graining in Fourier space (see the appendix), where all matrices are 
now scalars dependent on wavenumbers denoted by k or k, and frequencies denoted by w. We may formally rewrite 
the iterative relation for L in Fourier space as, 

L^G(^) = La{-, 2^^^^ ^^^2 ' 2 ^ ^^ /^*^ + ^^^ 2 ' 2 ^ ^^^ ■ ^^^^ 

Each successive coarse graining procedure gives us a new operator which weighs information from two different points 
of Fourier space, corresponding to wave modes of different length scales, and puts them into a new point. Even though 
the original linear operator contains only differentiation in space, the new linear operator after one step of CG has 
a time differentiation component as well. For uj — Q^ we can prove analytically (and verify numerically) that the 
operator reaches a fixed point. 



i(fc)=(^sin^^/(l-^sin^^). (16) 

This is the perfect operator for —d^ in one dimension. One might hope that this operator can be recombined with dt 
and used in the dynamic equation to give a perfect dynamics. It turns out that this is in general mcorrect. The reason 
is that the iterative relation from the path-integral calculation is a dynamic iterative relation with time derivative in 
it. When one sets w = 0, physically it translates into the assumption that small scale degrees of freedom are enslaved 
by the large scale dynamics. The small scale degrees of freedom instantaneously adjust to the large scale ones which 
are kept after each magnifying transformation. This is not physical. 

Since we are only magnifying in space, the time differentiation is diagonal in this phase space. We have the trivial 
relations, {dt)A = {dt)B = dt and {dt)c — {dt)D ~ 0. We define the full space-time evolution operator, 

L^ = dt + L such that L^cj) — V i (17) 

and the action operator, 

H = iZpL^ such that / Dry exp{-— ^ ryp?/} (5(?7 - L^</>) = exp{-— ^ ^iJ^} , (18) 

and express the iterative relation in terms of L^ and H . This leads to a simple form for the full iterative relation (see 
appendix) , 

where the constant factor c is 1 for CG and 2 for US. The second iterative relation physically means that the coarse 
grained version of the two point function of the true dynamics is preserved, if the coarse grained variable is governed 
by the operator L^^ with a non-Markovian noise source p. The above iterative relations are readily generalized when 
magnifications are carried out along the time direction. 

We now wish to determine the fixed point solutions of the operators L^ and H under their iterative relations. It 
can be shown that the operators approach their fixed points exponentially fast as a function of the number of iterative 
steps, irrespective of their detailed form at the microscopic scale. The fixed point solutions are given below while the 
exponential approach is illustrated in Figure w. 

We begin the the simpler case of US. Starting from a zeroth order operator of the form i(^.o — iuj + ^ sin^ (I), 
appropriate for a description at the microscopic scale e, after repeated US transformations we arrive at the operator 
suitable for the length scale Ax„ = 2"e. If the general form of the US operator after n iterations is written as 

it is closed under iteration, given starting values ao — (3o — 1- The iteration relations are 

c^n+i = a„ ( 1 + /3„ — ^ j and /3„+i = /3„ ( 1 + a„ — ^ j . (21) 

These have a fixed point solution 



Oiri 



i + i(*e„) + ^(ze„)2 + ... = j2mTTvM^nY = ;^sinh(^Ae;) 



Pn = I + M^On) + ^{tOn? + ■ ■ ■ = E(27w(»0")^ = -|;^ (cOsh( V^O;:) - 1) 

where 9„ = wAx^. 
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FIG. 1. RG flow of the dynamics operator. In this case the starting point is a microscopic Laplacian operator of the form 
I/o,tj = icj + ^fc^. The functional form of the n"^ iterate of L^j is, {Luj,n) 



2" i.^;3„ + _!-2-/„(fe)- 



For the iterative process starting from, for instance, the microscopic action operator iJg = cj^ + [^-sin (|)1 , we 
take 



if„ — 



a„ + e„sin^(|) 



62c.2 + (d„c. + ^sin2(|))2' 
where d„ and 6„ are the real and imaginary parts of /3 in L^.u- The iteration relations for a„ and e„ are 

fln+i = On + 2a„d„— ^ + (2a„ + e„)(6^ + c'i;) f -^ 



(22) 



(23) 



The fixed point solutions are, setting 6'„ = ■\/8„/2, 



-1 , o^ . 

"" ^ -•- + 360 ^ 



= E 74?f3V [4'+'] - ie„ = 44r(sinh(20„) - sin(20„)) - ie„ 






(4i+3)! 



3 + 7! + ^2.^ (4i+3)! ^^^ -^/ \ 



-^ 360 
6 ^ 



e 



E(4^[2(-l)^ 

E (4i+4)! i2(-l) 



i+ll 



^(sinh(6'„) cos(6'„) - cosh(0„) sin(6'„)) 
^sinh(6'„)sin(6'„) 
^(cosh(6'„)cos(6'„) - 1). 



(24) 



We can now move on to the CG case. Here we parameterize the operators as 



1 ^x„ 

^^.n — In -. r - 



4 ia;/3„ + s|jsin^(|)' 

rr-i ^ f (^.2 , an + e„sin^(|) 



(25) 



It is easy to see that CG shares the same /3, b and d parameters with US. The iteration relations for the other 
parameters are different. However, one can obtain a relation between a and a , namely, 



Using this relation we find 



^CG „US fj 1 ( „US ^US\ /0R\ 






Therefore, the fixed point solution for a'~^'^ and 7*-^*^ can be written in terms of that of a^^ , while the rest of the 
parameters have fixed points 

an =1- ■& + ••• =E(4ife)T[2''+^(-R-^-l)]-ie„ = (cosh(0„) - cos(0„))2Z„ 

e„ = -1 + ^ + • • • = E {^[8 ~ 2*'+7] ^ ^(^^^^^ _ ^„) 

/" = * - ^ + • • • = E (4i^[16(-l)i = -^(Z„ - 1) , 

where Z„ = 2|-[cosh(0„) sin(6'„) + sinh(0„) cos(0„)], and RI denote the average of real and complex parts of (^^)'*"+^. 

B. Perfect Action Operator in Space-Time and Stochastic CG Scheme 

So far, we have only coarse grained the spatial degree of freedom and obtained the corresponding perfect operators. 
In order to move on to numerical calculations on a lattice, we also need to coarse grain the time degree of freedom. 

We focus on the perfect action operator H = L^ pL^^ which is used later in the space-time Monte Carlo calculations. 
Here we derive the fixed point solution of H. We give a nearly closed form solution for H{k,u!) and show that this 
operator gives a perfect dispersion relation as measured from the time displaced two point function. The stochastic 
coarse graining scheme is introduced, which modifies H to give us an operator with reduced range of interaction. 

The iterative relation we developed previously does not hinge on whether CG was carried out on the space or time 
axis. Therefore, we can use it to CG in the time direction as well. Either one can start from a continuous description 
and alternately CG in space and in time, or one can directly use the perfect operator we developed previously and 
only CG from continuous time. Now, there is another dimensionless parameter, namely the ratio of the time scale over 
the characteristic time appropriate for a chosen length scale. For the diffusion equation, it is Ai/Ax^. We already 
see the manifestation of this parameter in the perfect operator derived earlier, where only the combination of the 
form wAx^ enters the expressions. Therefore, there are two restrictions on how we apply the two schemes. In the 
first case, we should CG twice in time direction for each CG operation in the spatial direction, maintaining the value 
of the ratio At/Ax'^ throughout the process. This means, for any reasonable values of At/ Ax^ at the macroscopic 
side, we need to start with a small Ax and a very small At. In the second case, we will not be able to maintain the 
ratio of At/Ax^. Therefore, the fixed point operator should be identified by iterating backwards. This means, we 
repeat the iterative process many times starting from various values of Ai„ = Ai/2" and iterate n steps. The fixed 
point is identified as the operator which is (within tolerance) not changed whether we start from At„ or Ai„+i. This 
method was used in the previous section to calculate the fixed point operator form for H when the time frequency lo 
was nonzero. This reversed iteration scheme is more powerful, since it can be generalized to other cases where there 
are other dimensionless parameters, such as in the case of massive fields. 

The fixed point solution of a d-dimensional operator under the CG iterative relation can be found using the 
techniques that have been described in this paper. An alternative method, the so-called "blocking from continuum" 
can also be used. In either case one finds 120|-|22^ 



where 0(p) is the continuum spectrum of the operator and 1 is a vector whose elements are of all possible integer 
values. 

In the above equation, an extra constant term with a parameter kq is introduced. This term is important for 
obtaining a localized perfect operator fit for numerical simulations P2|. To get this term, we modify the CG procedure 



to be a stochastic CG operation, also called soft CG, instead of hard CG, where an artificial noise term is introduced 
into the CG variable, 

0^ = + 1^, (29) 

where < v >= and (j/(i) i^(i')) = ^ ^y '^i.i'- Taking kq -^ oo, the hard CG case is recovered. 
Now consider the diffusion equation for a massive field, 

dt(t) = dl(j>-m(j> + r]. (30) 

The continuum spectrum of _ff is, 

i?=(|^)' + [(^)'+m]2 where co,k e {~n,7r) . (31) 

Defining the notation xi — x + 27r/, we have, 

1 i_^ 1 4sin^(fc/2)4sin^(a,/2) 1 



Aa;* ^ (A:;2 + /i)2 + r^ u;f, kf uif, 3 r^ k ' 

where we defined parameters /i = mAx^ and r = Ax^/Ai. To conform with notations used in quantum field theories, 
we have defined k = Kq At'^/3. 

The double summation is cumbersome to evaluate numerically due to its power decaying behavior. By rewriting the 
factor {kf [{kf + /i)^ + r'^ujf,)]}~^ as a difference of two terms we can re-express the above formula as a sum of a closed 
formed expression and an exponentially decaying expression. To do so, it is convenient to introduce w* = {kf + /i)/r 
and the function 

y^ 4sin^(fc/2) _ 1 (sinh^)(l-cosfc) 

^ '^'~ ^ kf{kf+fi) /i^ 7A^((cosh7M)-cosfc)^^' 

Then, after some simple algebraic manipulation, one finds 

-i^g-^ = -9,G-rsin^(^)9gG + 2rsin^(^)V f,tt^'/^) ^7' 7^"^" + ^ . (33) 

Aa;4 ^ ^ 2 ^ '" ^ 2 ^ ^^ fc,2(fcf + ^)3 coshtjf - cost^ 3 r2 k ^ ' 

Now what remains of the summation is much easier to evaluate due to its exponential decaying behavior. 

From the above equation, we can obtain the dispersion relation implied by such an operator. The two point function 
for a free field is S{k,uj) = H~^{k,uj). Taking the discrete Fourier transform of equation (p3) back to real time gives 
the static equal time structure factor, 



5,., t.O)^j: '-^J^i^! - 1 + .-"■•) + 3^ , (34) 

and the time displaced two point function 

-^ kf 2(c.;)3 ■ (35) 

All dynamic modes are present, each with the correct decaying behavior and with a prefactor (enclosed in curly 
bracket) due to coarse graining in space as well as in the time direction. In principle, the decay rate should be 
measured in the long time limit where all modes outside the first Brillouin zone are negligible. However, for all 
practical purposes, the / ^ modes are negligible (or more precisely, the next significant mode not degenerate with 
I = 0) even for short times. For example, for k = 7r/2, fi = 0, the amplitude of the next most significant mode 
{I = —1) is only 1.5 x 10"'' of that of the I = mode. Therefore, we can use the i > 1 values of the time displaced 
two point function to evaluate the perfect dispersion relation for all the wave modes with wavenumber within the first 
Brillouin zone. 



From H~^{k,LL!), we obtain the perfect operator coefRcients H{r,t) in real space and time. Notice that "the fixed 
point of an operator" actually means the fixed point of the dimensionless operator. Consequently, operator coefficients 
for the perfect action operator are actually those of H Ax'^. For practical reasons, we need to adjust the parameter 
K for optimal locality. In one dimension, k « 2 and 6 are the best values for d^ and d^ respectively. Therefore, we 
need to find a compromise. The best scheme is to choose k = 2 such that the most significant couplings lie within a 
rectangle area elongated along the x direction. This way, the total number of significant couplings is minimized. 

The leading order coefficients of H for k ~ 2 and zero mass are tabulated in Table and are shown in Figure 0. 

TABLE I. Sample coefficients of the perfect action operator for tlie diffusion equation, k = 2, fi = and Ax'^/At = 1. 



{t,x) 


H 


(t,x) 


H 


{t,x) 


H 


(i,x) 


H 


(0,0) 


3.90458 


(0,1) 


-1.02978 


(0,2) 


-0.0421266 


(0,3) 


0.098042 


(0,4) 


0.0291451 


(0,5) 


-0.00317113 


(0,6) 


-0.00407848 


(0,7) 


-7.35334 X 10"* 


(1,0) 


-0.464966 


(1,1) 


-0.278677 


(1,2) 


-0.0339122 


(1,3) 


0.0328692 


(1,4) 


0.0148371 


(1,5) 


-2.8286 X 10"* 


(1,6) 


-0.00204057 


(1,7) 


-5.22915 X 10"* 


(2,0) 


-5.99324 X 10~* 


(2,1) 


-8.25418 X 10"* 


(2,2) 


-3.68371 X 10"* 


(2,3) 


6.91673 X 10"* 


(2,4) 


8.04927 X 10"* 


(2,5) 


2.04135 X 10"* 


(2,6) 


-1.27393 X 10"* 


(2,7) 


-9.25806 X 10"" 
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FIG. 2. Surface plot of amplitude of perfect action operator coefficients for the diffusion equation. The coefficients exponen- 
tially decay away from the origin. The decay speed is slow along the x direction, k = 2, m = and Ai = Aa;^. 

IV. NUMERICAL SIMULATION USING THE PERFECT OPERATOR 

In this section, we discuss the application of the perfect linear operator in numerical simulations of Langcvin 
dynamics. We show that the perfect operator should be decomposed into an Up operator and a Down operator in 
order to obtain a correct equation with a finite number of high order time derivatives. Without this decomposition, the 
truncation of the perfect operator is highly non-trivial, if not impossible. For Langcvin dynamics, the dynamics of the 
non-Markovian noise is difficult to obtain because it requires taking the square root of the noise correlation function. 
Various numerical simulations were carried out using the truncated perfect operator and other approximations, to 
illustrate the advantage of using coarse grained variables as opposed to uniformly sampled variables in numerical 
simulations. These, together with the limitations of this approach, are also discussed. 

A. Perfect Operators in Langevin Dynamics 

Here we derive the perfect operators U and g appropriate for Langevin dynamics. In the previous section, we 
obtained the iterative relation for both L^ and H . From the latter, the correlation function for the non-Markovian 
noise is obtained. The discretized system follows the dynamics described by the PDE, 



i„ 



■n, 



(36) 



where rj satisfies {ri{i,j)ri{i',j')) = 'SyP ^(* ^ i\j ^ ]')■ This formally simple equation is different from the usual 
Langevin dynamics in two respects: the non-Markovian nature of noise and the presence of (in principle) infinite 
orders of time derivatives in both L^^ and p^^ . 

The non-Markovian nature of the noise means that there is dynamics in the noise variable. This is not surprising. 
In the path-integral calculation, each CG step results in formally discarding small scale degrees of freedom. But in 
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fact, the small scale degrees of freedom are not entirely discarded. Since the small scale dynamics is affected by the 
noise source as well as the system dynamics at the coarse grained level, when the small scale degrees of freedom are 
integrated out at each CG step, part of the small scale dynamics is preserved by modifying the dynamics at the larger 
length scale and by injecting dynamics into the noise. This is basically a feedback effect. 

Due to the non-Markovian nature of the noise, we need to write down the dynamics followed by the noise, 

p'^'v = Vo , (37) 

where 770 is a white noise satisfying {'r]o{i,j)riQ{i',j')) — -^6i,i' Sjj'. The matrix p^" is the square root of p in the 
sense that the product of p^'^ and its Hermitian conjugate gives p. For instance, in Fourier space, Vw^ + fc^ = iuo + k"^. 
There are in principle infinite orders of time derivatives in p^'^, just as in L^^. 

Naively, L^ can be obtained as a series expansion in dt which is then truncated to certain order. This turns out 
not to be the correct approach. Rather, we need to decompose the operator L^ in the form of a numerator {U) over 
a denominator (-D), 

L^ = D-^U, (38) 

where we write the denominator as an inverse operator. The distinction between the numerator and denominator is 
easily seen in the fixed point operator. We can eliminate the inverse operator by applying D on both sides of equation 
(pq). Redefining the noise as ^ = Dj] and denoting its correlation function g~^, we have, 

U-'g-'tf^-'^L-'p-'L'^-K (39) 

The operators U and g are therefore equivalent to the older pair of L and p in the evolution of the discretized system. 
Equation (Bq) may be rewritten as 

f7^0 = e with g^/^^ = VQ- (40) 

Using (p5|), and in the notation of the last section, the perfect operators for the diffusion equation under the CG 
scheme are 



U^ = zw/3 + ^-2 sin^ - , 

ei/^ = {a + 7sin^^+7/3f}/{a + esin^^ + /6^(|)^ + /(d|+sin2^)^}^. (41) 

For analytical tractability, we used the closed- form solutions of the operators available for discrete space but continuous 
time. 

Unlike L, the new evolution operator U can be expressed in a clear and simple series expansion. The spatial part 
is simply the central difference operator and the time part is a sum of all orders of time derivatives with constant and 
fast decaying coefficients (see equation (|2^)). 

The operator p^ " has a very complicated form. It has many high order space and time derivatives, which in general 
are coupled. Series expansion and truncation are necessary. To the first order in Ax^, we have for CG, 

g,i/2^ l_i sm^{k/2) + ciujAx'^ = I + X^ dl - rdt , (42) 

6 

where c — 1/6 — l/\/720 = 0.129, A = Aa;/\/24 and r = cAx^. Therefore, the noise source is largely a white noise. 
It has a correlation length of the order A and a relaxation time of the order r. When the form of the operator is 
obtained and truncated to a specified order, one can evolve the system according to equation (HQ). 

Often periodic boundaries are used in the spatial dimensions. Therefore high order spatial derivatives do not pose 
a problem on a lattice. Higher order time derivatives, however, require a corresponding number of initial conditions. 
This might pose a problem, especially for the non-Markovian noise. If one is interested in equilibrium properties of 
the system, the initial transient stage is not important. An initial condition with all derivatives zero is fine. When 
one wants to study the initial transient stage corresponding to a certain microscopic initial condition, one can evolve 
the system using a fine mesh for n steps under a conventional numerical scheme, where n is the highest order time 
derivatives. For each step, one can coarse grain the microscopic configuration to the desired CG level, insert the CG 
version of cj) into equation (^^j, and the noise in the transient stage is obtained. This way, initial time derivatives for 
both coarse grained and ^ can be computed. 

The calculation of the space-time discretized g^''^ can be quite involved [ p3[ . Since our main interest is in calculating 
equilibrium properties of dynamic systems, we can take an alternative route, namely Monte Carlo simulation, as 
discussed later. In this case, the perfect action operator H is all we need. 

12 



B. An Example of Using Perfect Operator in Langevin Dynamics 

In this section, we present an application of the operator U to the deterministic dynamics of the coarse grained 
variable governed by the diffusion equation. The (truncated) perfect operator U gives superior results for the evolution 
of the configuration. The relative advantage of using the CG variable vs the US variable is also touched upon and 
will be studied more closely in the ensuing section. 

For simplicity, we truncate the series expansion of (3 to the first order to obtain an operator U with a second order 
time derivative. Direct truncation is not appropriate when setting higher order terms to zero, since we should adjust 
the remaining coefficients. Instead, we use the operator at the first level of CG, starting from a central difference 
operator. The coefficients for time derivatives higher than the second order are identically zero. We have, 

Ax"^ 1 

-— a,^</., + M, + -^^(2 0. - 0,+i - (/..-i) = 0. (43) 

lb Aa;^ 

Suppose the system is periodic with length L. The initial condition has modes down to length scale e = L/M with 
AI being an integer, namely, 

M/2 

0(m,t = O)= Y. e^'"""'^'^'/'fc , (44) 

where (/)fc is the amplitude of the kth wave mode. We know analytically the exact solution: by coarse graining the 
exact solution to a length scale Aa; = L/N = pe, we have, 

N/2 p/2 ■ r A /T\ ^/2 

q=~N/2 i=-p/2 ^ \ \'i ' J I J q=-N/2 

where <f>q,o{t) is the exact wave mode for the CG variable. This equation gives us both (j){n,t = 0) and dt(j){n,t = 0). 
Now let us ask: what result would equation (|3|) yield on a lattice with grid size Ax, given the CG initial conditions? 

We have 0g(i) = I]f^ip/2 Cq^,{t) (j)q+iN, where 

psm(7r(g + I A')e/L) Acj L'^ 

(jj- — 3^ sin i-^), and Acj = ■^^cos(-^). For comparison, the corresponding result from conventional numerical 
analysis (NA), which is the same as just keeping the first-order time derivative in U , is 



.NAu^ r 4sin"(7rg/iV)^ 

Aa;2 



C^At)^eM- -pj-> t}, (47) 



where the time evolution does not depend on i. The solution for modes within the first Brillouin zone, i.e. i = 0, is 
greatly improved as shown in Figure 0, where we have plotted the time evolution of the coefficient Cq^i=o (i) (without 
the prefactor due to CG) for selected q values. For small q, the Aw dependent part in equation (|6|) is not important. 
A Taylor expansion clearly shows that cj_ is closer to the true decay rate than the NA result. 



13 




Exact. q=7t/3 

o o NA. q=3i/3 
PO. q=7i/3 
+ + Exact. q=27i/3 
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FIG. 3. Decay of wave modes in the first Brillouin zone using PO and NA equations vs the exact result. The decay rates for 
the PO scheme are closer to the exact ones than the NA results. The coefficient is Cq,i=o(t) without the CG prefactor. Sample 
wavenumbers are q — n/3 and 2tt/3. 
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FIG. 4. Decay of wave modes in the second (i = —1) and third (i — 1) Brillouin zones using PO and NA equations vs the 
exact result. For the PO and NA schemes, wave modes do not decay as fast as for the exact result. PO is better for modes in 
the second Brillouin zone than NA and is also advantageous for late times for modes in other Brillouin zones. PO results are 
very close to the exact one at short times as indicated by the dip in the plotted curve. The coefficient is Cq,i^o{t) without the 
CG prefactor, which reduces the importance of modes outside the first Brillouin zone. Sample wavenumber is g = 7r/2. Notice 
that the NA result does not depend on i. 

For very short times, the dynamics of all modes are correctly prescribed, even for i 7^ 0. This manifests itself as 
a dip in the short time region of Figure ^ (subject to resolution in the time axis, the dip for the i = 1 mode is not 
discernible in the figure). For finite time, modes outside the first Brillouin zone decay quickly in the true dynamics 
(Figure Q). In the PO result, the decay rate is dependent on q, therefore these modes do not decay as fast as they 
should do. But since the PO result also contains information on i, for modes in the second Brillouin zone, the resulting 
dynamics are still closer to the true one than the NA result. This is because we used the PO operator U of one level 
CG. For higher wavenumbers, due to the i dependent term, there is an anomalous (negative) amplification of wave 
modes at the initial transient stage which disappears later. Therefore, the power spectrum of the configuration should 
die off quickly for modes with length scale much less than Ax. In other words, we should not over-coarse-grain. It 
follows that as we keep more and more terms in the perturbative series for U^j, the PO result will be close to the 
exact one for higher and higher wave modes. 

The prefactor 

sin(7rgAa;/-L) 
sin(7r(g + i 7V)e/L) 

modifies the contribution of each wave mode to the solution. This comes from using the CG variable in the PDE and 
is very important in reducing errors that arise from using the discretized PDE. For instance, although modes with 
g « decay very slowly in the PO result, their prefactor is close to zero for i 7^ 0, while they very quickly decay to 
zero in the true dynamics. 

Notice that US and CG share the same tj . In the US scheme, there is no prefactor. Modes with g w and i ^ Q 
do not decay. If we use the same equation as above, the prefactor for t > Ax^ is, 



y 



C,.o,.(i) «!-(-)'■ (48) 



For large i, it overstates the contribution of the mode to the solution and is worse than NA. This imposes a stricter 
constraint than for the CG scheme on the power spectrum of the configuration, and is the reason why CG is a better 
scheme. This has been tested numerically on several model dynamics |l[. 

V. SPACE-TIME MONTE CARLO SIMULATION 

The path-integral formulation easily leads to a space-time Monte Carlo simulation. We discuss issues related to 
truncating the perfect operator such that it has a finite range of interaction. Numerical simulations are carried out on 
the linear diffusion equation to test computational efficiency of using the perfect operator, and on Model A dynamics 
to test the merit of direct application of the perfect linear operator to nonlinear dynamics. 

In quantum field theories, many problems are formulated in terms of path integrals. Numerical simulations usually 
employ the Monte Carlo method, where due to space-time symmetry, time is simply treated as one of the dimensions 
in a d -|- I-dimensional lattice. In statistical physics, when dynamics is involved, evolving a Langevin equation is 
the norm. A typical form of the equation contains a first-order time derivative, a diffusion term and some nonlinear 
interaction. Time and space are not symmetric. However, numerical simulation of a Langevin equation is not the 
only choice for studying dynamics. We can also perform Monte Carlo simulation on a space-time lattice [gj], similar 
to the approach adopted in quantum field theories. The basis for such a calculation is the path-integral formulation. 
Starting from equation (0) , and performing a trivial integration over the noise to eliminate the delta function, we have 

P=fD<p exp{-^ ^[(^0 + fm^ - ^d^f]} . (49) 

The cross term linear in dt results in a boundary term and does not influence directly the calculation of P. Ignoring 
the Jacobian contribution, we are left with a positive definite functional. We call the term in the exponent the 'action' 
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for obvious reasons. For linear operators, we know how to coarse grain the above expression, 
noise in the above equation, we have 

P= fi^^expl-^^^iJ^}, 



Integrating out the 



(50) 



where H is the fixed point operator of the action operator. Working with this path-integral formulation, we do not 
have to worry about taking the square root of the noise correlation matrix as we would with the Langevin equation. 
In the following, we will look at a specific example of the linear theory, namely the dynamics of a system described 
by the diffusion equation. 



dt (p ~ dl (f) — iTKp + rj , 



(51) 



where tti is a constant which we will call the mass and t] is white noise with strength fJ. We have chosen a unit 
diffusion constant. In the space-time Monte Carlo probability we use the 1 + 1-dimensional perfect operator for 
—dl + (—9^ -f m)^ developed in the previous sections. Then we look at the application of the perfect linear operator 
to the nonlinear Model A dynamics. 

A. Truncated Perfect Operator 

The perfect operator needs to be truncated to finite range to be used in numerical simulations. Although the 
introduction of stochastic CG reduced the interaction range of the perfect action operator, the operator coefficients 
do not terminate in a finite range. Furthermore, they decay slowly along the x direction, where the coefficient of the 
10th neighbor |p5[ still has an amplitude of around 1 x 10~^. This makes truncation of the perfect operator more 
problematic than in quantum field theories, where keeping next nearest neighbors is already very good ||22| . 

One criterion for truncation is that the magnitude of the discarded coefficients have to be small. But there are other 
considerations as well |2^,^. One would like the operator to satisfy certain constraints that stipulate the correct 
behavior of the operator in the continuum limit. These constraints are in the form of sum rules p2]. For the diffusion 
action above, the constraints in the continuous limit are. 






2 Z^i 1^3 -"».: 



-2fi 



(52) 



^,3 



= 1, 



where, as defined previously, /x = m Ax^ and r — Ax'^/At. Naively, one might expect that one way of proceeding 
would be to truncate the perfect operator to a finite and manageable range, and then to enforce these constraints 
to improve the directly truncated operator. In reality, these sum rules are not satisfied even for the perfect operator 
for finite Ax and finite n. The error is of higher order in Aa: and inversely proportional to k. On the one hand, the 
continuous limit constraint conditions can be recovered. On the other hand, for finite Ax, the constraints no longer 
hold unless an operator with a long interaction range is used. The average constraint error is about 0.1% if one keeps 
up to ~ 20 and '^ 3 neighbors in the x and t directions respectively. 



TABLE II. Coefficients of naturally truncated 11x3 perfect action operator for diffusion equation, k = 2, fi = mAx^ = 0.25, 
Ax^/At = 1. 



it,x) 


H 


(t,x) 


H 


{t,x) 


H 


{t,x) 


H 


(0,0) 


4.00869 


(0,1) 


-1.00198 


(0,2) 


-0.0819891 


(0,3) 


0.0724608 


(0,4) 


0.0270167 


(0,5) 


3.546940 X 10"* 


(0,6) 


-0.002564 


(0,7) 


-7.593954 x 10"* 


(0,8) 


4.298800 X IQ-" 


(0,9) 


8.703907 X 10"'' 


(0, 10) 


1.817881 X 10"" 


(1,0) 


-0.430984 


(1,1) 


-0.265854 


(1,2) 


-0.046095 


(1,3) 


0.021651 


(1,4) 


0.012848 


(1,5) 


0.001211 


(1,6) 


-0.001157 


(1,7) 


-4.724628 x 10"* 


(1,8) 


-6.542402 X 10"" 


(1,9) 


4.851149 X 10"'' 


(1, 10) 


1.406475 X 10"" 


(2,0) 


3.220461 X 10"* 


(2,1) 


-2.491026 X 10"* 


(2,2) 


-5.120127 X 10"" 


(2,3) 


1.844663 X 10"* 


(2,4) 


5.223556 x 10"* 


(2,5) 


2.386099 X 10"* 


(2,6) 


-2.442266 x 10"'' 


(2,7) 


-5.894465 x 10"" 


(2,8) 


-1.676697 X 10"" 


(2,9) 


3.755505 X 10"" 


(2, 10) 


4.507152 X 10"'' 
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FIG. 5. Decay rate of wave modes for diffusion equation, k = 2, fi — 0.25 and Aa;^/Af — 1. Perfect operator decay rates are 
obtained using the first two t 7^ nodes (equation (3E|)). 



An alternative approach [g6|,^ is to compute the perfect operator on a smaller lattice and then use this 'naturally' 
truncated perfect operator. In this way, the constraint is taken care of in the continuum limit. If we use the operator 
on a lattice of the same size, the operator gives a perfect dispersion relation. However, when it is used on a larger 
lattice, it is no longer perfect, as can be seen from the inexact dispersion relation for high wave number modes, which 
are those most affected by truncation. The reason lies in the high decay rate associated with a k^ dispersion relation. 
For k = TT, the ratio between successive S{k,t) values is about 2 x 10^. Thus, to maintain exponentially decaying 
scaling over three nodes, we need a relative accuracy of 10~*. Taking into consideration the importance of keeping 
enough neighbors and the computational efficiency, an operator with up to 10th and 2nd neighbors in the x and t 
directions is chosen as the operator for most of the subsequent computer simulations. An operator with 9th and 2nd 
neighbors in the x and t directions is also used for some of the simulations. There is no discernible difference between 
this operator and the 11x3 one. 

The operator coefficients are displayed in Table || for ^ = 0.25 and At = Ax'^. For the above operator with 
/i = 0.25, a 3 node scaling regime is maintained for 60% of the k mode and a 2 node scaling regime for about 94% of 
k mode. For a larger operator of size 20 x 7, we would have a 3 node scaling regime for about 90% of the k mode. 
The decay rates for different operators are compared in Figure g. 

The rapid decay rate of high wavenumber modes is what distinguishes the perfect operator for the diffusion equation 
from the 1 + 1-dimensional Laplacian operator used in high energy physics. In the latter case, the ratio between 
successive S{k, t) values is at the more benign level of about 0.04. The exponential decaying range spans more values 
of time displacement. It is easier, therefore, to read off the dispersion relation all the way to the edge of the Brillouin 
zone. It is also more stable with respect to small changes in coefficients of the operator. 
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B. Numerical Simulation of Diffusion Equation 

We carried out space-time Monte Carlo simulations to test the efficacy of the perfect operator developed in the 
previous section. Suppose we are interested in the diffusion dynamics of the system described by equation (^l|) and 
would like to calculate its space-time correlation function. Let the system be of length L = 16, with the spatial scale 
of interest I — 1. In the path- integral formulation, the time span of the system is T = 8. Both space and time 
directions have periodic boundary conditions. The Metropolis algorithm is employed [ pSj . 

Three simulation runs are presented. One simulation uses a perfect operator with range of interaction up to 10th 
and 2nd neighbor in the x and t directions respectively. A lattice of Nx x A^f = 32 x 32 was used, corresponding 
to Ax — 0.5. The other two simulations were carried out on 32 x 32 and 64 x 128 lattices using the conventional 
central difference operator. In each case, the time direction grid size is Ai — Ax^. In each simulation, Nrun number 
of independent runs were conducted to obtain statistics of measurements, each run with A^ = 5 x 10^ MC steps (one 
sweep of the system) and one measurement per 8 steps. Nrun = 6 and 7 for the 32 x 32 and 64 x 128 lattices 
respectively. The (fc ~ 0,^; ~ 0) modes have the largest standard error, which is crucially dependent on the lattice 
size. The typical percentage standard error of S{k,uj) for 32 x 32 lattice is about 1% and 2.5% for PO and NA 
operators, while that of 64 x 128 lattice is 6%. 




Exact 
o NA.Ax=0.5 
+ NA. A x=0.25 

PO.Ax=0.5 




FIG. 6. Cross sections of S{k,uj) for the diffusion equation, m = 1, L = 16, T = 8. Cross sections are at a; = (left) and 
fc = (right). The exact result is [(m + P)^ + uj^]~'^ . 

In Fourier space, cross sections of the space-time displaced two point function 5(fc,a;) are plotted in Figure y. We 
do not expect the perfect operator result to be exact because S{k, uj) is now a two point function of the CGed variable, 
not the continuous variable. But it turns out to be quite close to the exact result. The NA result for 32 x 32 lattice 
deviates further from the true value at the same (fc,w) value. For this plot, a constant offset of ^pj^^- is subtracted 
from S{k,uj) of the perfect operator runs to eliminate the contribution from the added noise in the stochastic CG 
transformation. 

Fourier transforming S{k,uj) to real time, we obtain the dispersion relation from S{k,t) ^ g-'^(fe)t_ Xo avoid static 
contributions in the t = mode, we choose the most significant t y^ points to calculate U!{k) — {log S{k, At) — 
logS{k,2At))/At. The results are shown in Figure M. The perfect operator gives a near 'perfect' dispersion relation 
for the length scale we are interested in (corresponding to wavenumber fc ~ tt), giving the correct zero k mode mass 
and correct k^ dependence. We can get a comparable result using a larger lattice with the NA operator, but with 
more computational effort. For large k modes, the amplitude of S{k, 2At) is of order 10^^ relative to that at i = 
and becomes unreliable given the simulation accuracy. The real value is used in the plot when S{k, 2Ai) is negative. 
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FIG. 7. Decay rate of wave modes for the diffusion equation, m = \^ L = 16, T = 8. Lattices yield At = Aa;'^. Length scale 



of interest corresponds to fc '^ tt. 
PO (a) are obtained using the t = 



Exact result is m + fe . PO results use the first two t 7^ nodes of S{k,t). NA results and 
= and t — At nodes. 



One might ask: why caU the operator perfect when it docs not reproduce the correct dispersion relation for 
wavenumbers beyond k = tt7 The answer is that it is not the operator that is not perfect but the simulation itself! 
The perfect operator gives the best result possible for physical quantities of interest given the error of the simulation. 
With more statistics, the dispersion relation from the perfect operator approaches the correct result for all modes 
with length scale larger than the grid size. The same is not true for the NA operator. For a discretization twice as 
fine, with increasing number of statistical samples, the dispersion relation for the NA operator approaches a limit 
that is different from the true solution, and is about 19% off at the edge of the Brillouin zone. 

The simulation error can be overcome when we choose smaller At relative to Ax^. As shown in Figure ^, the PO 
decay rates using At = i Aa;^ (corresponding to 32 x 64 lattice) closely follows the exact result and is more accurate 
than the measurement from NA. 
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FIG. 8. Decay rate of wave modes for the diffusion equation, m = 1, L = 16, T — %. Same as in Figure fTlexcept that lattices 
with At ^ Aa;'^ are used. 

One might as well choose operators according to the magnitude of the statistical error of a simulation. Given the 
usual error of 1% for S{k,uj), a smaller sized perfect operator could be used to improve efficiency of the simulation 
without compromising accuracy of the physical measurements. Even with a 11 x 3 PO as used in our simulation, the 
extra computational effort is not that huge. This operator requires 21x5= 105 points be used to calculate the action 
density at each grid point, whereas 7 points are used in conventional NA calculations. However, since most of the 
computation effort goes to generating random numbers (we used Numerical Recipe's ran2() subroutine p9[ as well as 
SPRNG modified lagged Fibonacci generator from NCSA JS^I), it turns out that the overhead from extra neighbors is 
not significant considering the improvement of results. If one uses a naturally truncated 5x2 PO, total CPU time for 
the sample calculation will be reduced by 58%. The decay rate rivals the result from NA with a lattice twice as large. 
In this case, however, we will not recover a perfect decay rate with better statistics due to the severe truncation. 

Our code is written in C++. On a SUN Ultra2200, the run times are shown in Table III for a test run on a 32 x 32 
lattice with 10000 MC steps. For the same number of steps and lattice size, the PO calculation takes about 4 times 
as much time as the NA calculation. Their standard errors for decay rates are roughly the same if the same nodes 
are used. However, the PO uses the second and third nodes to calculate decay rates. Therefore the resulting decay 
rates have standard errors about twice the size of that for NA. 



TABLE III. CPU time of simulations on diffusion equation using PO vs NA. 32 x 32 lattice. 10000 Monte Carlo steps. For 
the same number of statistical averages, the standard error of S{k, uj) for PO is about half of that for NA. 





Action Calculation 


Random Number Generation 


Total Time 


NA 


6.6s 


15.0s 


30.1s 


PO 


105.4s 


13.6s 


128.5s 
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The relevant quantity regarding the computational efficiency is the total computational effort (TCE) needed to 
reach a certain level of root mean square (RMS) error S. This is defined as 



TCE^cNNtNx, 
where the speed factor, c, is 4 and 1 for PO and NA respectively. The RMS error S is given by 

S' = Sl + 5l, 



(53) 



(54) 



where ^i is the bias and 62 is the standard error. In comparing the efficiencies of PO and NA, we focus on the wave 
mode with k = n. 

For the naturally truncated PO, Si ^ 0.01% and is negligible. For a 32 x 64 lattice, 64K MC steps are needed to 
reduce 62 to 1% for /c = tt. Hence TCE = 5.2 x 10^ 

For NA and a large lattice size, we have. 



a b k'^L^ 

oi ~ T7T + T77 7 where a — — — -ttt 

Nl N^ 12 (m + fc2) 



and b 



jm + efT^ 
24 



(55) 



For instance, with L = 16, T = 8, m — 1, k = tt, one has a = 191.2 and b — 315.1. The standard error $2 is inversely 
proportional to VjV and is a function of the lattice size. Increasing the lattice size increases 62- However, increasing 
Nt also has the effect of improving the result, since smaller At relaxes the constraint on the statistical accuracy of 
the first few nodes of S{k, t). Let us assume that 



62 = s^^^n:n,^/Vn, 

where a and /3 are constant parameters. The minimization of the total computational effort yields. 



(56) 



2_ 2 + 2a + 2^^ 

^ ^ 2a +1 Ui 

^2^ 2 + 2a + 2/3, b 



2/3+1 'Si' 

— = (1 + l±^±A) ll^fo^ (l^fP . ^2(A^o,^.,o,iVt,o) .2 
No 2 Nx,Q Ntn 5 



(57) 



where the optimal Si — S/ Jl + i/^,g and where S2{Nq, N^.o, Nt^) is the S2 value for a lattice size {Nxfi, Nt.o) and 

with A'o Monte Carlo steps. For instance, with the above a and b values and a = f] — 0, to reach a RMS error of 1%, 
one needs N^ = 257 and Nt = 330. Given that S2 « 1.2% for N^ = 40i^, N^^ = 128 and Nt^ = 256, we expect the 
optimal N = 86K. Therefore TCE = 7.3 x 10^. If we have a = 1 and (3 = instead, the optimal values are N^ = 190 
and Nt = 423 and N = 254K. Therefore TCE = 2.0 x lO^*'. There is a factor of 40 improvement (see Table |%. The 



advantage of PO will be more pronounced in higher dimensions. 

The values of a and (3 are difficult to obtain. The values a = 1 and /? = are good approximations for the relevant 
lattice sizes, namely N^ and Nt of order of or bigger than 200. Notice that a large lattice size is most detrimental to 
the standard error of the small k modes. 

TABLE IV. Total computational effort for PO vs NA. One requires that the root mean error of uj{k) be 5 = 1% for k = n. 
Parameters are a = 1 and /3 = (see equation (p9)). 
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Nt 


N {xW) 


TCE (xlO'*) 
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190 
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PO 
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FIG. 9. The standard error of the decay rate of wave modes for the diffusion equation for NA using different lattice sizes. 
m = 1, L = 16, T = 8. Standard errors are normalized to A^ = 10^ Monte Carlo steps. 

In summary, we find that the perfect linear operator gives us the perfect dynamics of the various wave modes, given 
the errors of a numerical simulation. For the same lattice size and number of Monte Carlo steps, the PO scheme 
(with the 11x3 operator) is about 4 times slower relative to the NA scheme, where generating random numbers takes 
about 50% of the total computation time in the latter case. However, the computational effort in order to reach the 
same root mean square error for PO is on the order of 1/40 of that for NA. This will be more pronounced in higher 
dimensions. Moreover, a more severe truncation of the perfect operator is possible, given the inherent accuracy of the 
simulation, further enhancing the efficiency of the PO scheme. 



C. Numerical Simulations on Model A Dynamics 



In this section we study the application of perfect linear operator to the time dependent Ginzburg-Landau equation 
for Model A dynamics, 



dt(t>^d'l(t>-m(f>~g4)^ + 'q. 



The corresponding path-integral formula is, 

P^li^^expl-^^i^o + ^i]}^ 



(58) 



(59) 



m)'^)(j) and Si ~ 2 g (p^ i^d^ + 'm)(j) + {g (fy^)"^ — SF'^^ ^^'^ contributions from the linear 



where 5*0 ~ 4>{—dl + (—9^ 

and nonlinear terms respectively. 

For systems with nonlinear interactions, an exact analytical expression for the perfect operator is not available. 
The difficulty lies in the fact that the form of the continuous action is not closed under the CG transformation. New 
interaction terms are generated in reaching the fixed point of the discrete description of the dynamics. In general 
there is an infinite number of interaction terms of diminishing importance. In order to proceed, we need to make some 
approximations. In conventional numerical analysis, the form of the continuous action is used, where the Laplacian 
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operator is replaced by the central difference operator and local self-interactions are left unchanged. In analogy, we 
use the perfect linear operator developed previously for 5*0, while leaving the nonlinear self interactions unchanged. 
We bundle the rmp term in with the g 4>^ term in the to < regime to reduce the standard error of the numerical 
simulation. Intuitively this is a reasonable thing to do since |(/)| develops a non-zero amplitude and the contribution 
to the dynamics of (f) from these two terms largely cancel each other. We used the conventional central difference 
operator for the operator —d^ + to in 5*1. 

There are two regimes: to, > where the nonlinear term amounts to a renormalization of the mass, and tti < 
where a nontrivial ground state develops with a magnitude ±^Jm/ g. 

1. The m > Regime 

We simulated the dynamics of a system of physical lengths L = 16, T = 8 and parameters m — g = ^ — Ion lattices 
of different sizes. Mass dependent perfect linear operators are used. The Fourier transformed space-time correlation 
functions S{k, oj) are measured and averaged over several runs. Most simulations consist of Nrun — 9 runs, each with 
TV = 3 X 10^ Monte Carlo steps. Measurements are done every 8 Monte Carlo steps. For the NA result with 64 x 256 
lattice, 8 runs are used. Fourier transforming S{k,uj) to real time, we obtain S{k,t), where S{k,t — 0) is the static 
structure factor and the mode decay rates can be read off from the time dependence of S{k,t). The length scales of 
interest are those larger than Ax = 1. As in the case of the diffusion equation, the standard error of the PO result is 
half of that for NA with the same number of statistical averages. 

Mode decay rates obtained from the PO scheme for k away from the origin are greatly improved over its NA 
counterparts, as shown in Figure Kn For Ax = 0.5, At — 0.25, if we had used the second and third nodes of S{k, t), 
the decay rates for the second half of the Brillouin zone would not be reliable, reflecting the inherent numerical error 
(roughly 1%) of the simulation. This is as in the free field case discussed at the end of the previous section. For the 
plots, we used the t = and t = At nodes instead. It is no longer perfect, but it is within the numerical error of the 
simulation and gives improved results as compared with NA. When we choose At — 0.125, the error of the simulation 
is no longer a limiting factor and the decay rates over the whole Brillouin zone are recovered using PO. With a smaller 
At/ Ax^ ratio, the time direction becomes more continuous and the decay rate values are improved for all schemes as 
expected. 
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FIG. 10. Decay rates of wave modes for the Ginzburg-Landau equation, m = 1, L = 16, T = i 



For TO > 0, the ground state of the order parameter has an expectation value of zero. The nonhnear self interaction 
term in equation ( pq ) has the main effect of renormalizing the mass to a new effective mass m^s = m + g{(fP'). 
In mean- field theory, the expectation value of (jP' is expressed as a function of m^^s, which is then self-consistently 
determined by the relation 



rries/m = 1 



1 






The renormalized mass is easily seen to be larger than the bare one. 



(60) 
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FIG. 11. The influence of Ai/Aa;^ on the decay rate of small k wave modes for the Ginzburg-Landau equation. Smaller At 
gives improved result for NA. The effective mass for PO, however, approaches a limit less than the mean-fleld result, m — 1, 
L = 16, r = 8. 

From the decay rate of wave modes with fc « 0, we can read off the value of the renormahzed effective mass. The 
mean field value of the effective mass is nics = 1-2258 for the chosen parameters. For the NA scheme, the renormalized 
mass is less than the bare mass when the grid size along the time direction is chosen to be At = Ax^. Reducing the 
grid sizes while retaining the ratio At/Ax'^ leads to reduced effective mass values, away from the correct result. For 
the 64 X 128 lattice, we have nieS ~ 0.26. Unlike in quantum field theories, time and space are not symmetric in the 
dynamics we are considering. This translates into a freedom of choice of grid sizes At and Ax. Physical considerations 
lead us to the natural choice of At = c • Ax^ where z is the dynamic exponent and c is a constant factor. Outside 
the critical regime, the diffusion term dominates the dynamics and z equals the mean field value of 2. We expect the 
constant factor c to be dependent on the nature of the nonlinear interaction and to be different from 1. When we over 
coarse grain in the time direction relative to the space direction, the (relatively) finite size of At introduces error into 
the simulation results. We found that a Ax^/At ratio value of 2 to 4 is needed to reduce this error (see Figure O). 
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FIG. 12. The static structure factor S{k, t = 0) for the Ginzburg-Landau equation. 



1, L = 16, r = 8. 



For the PO scheme, the effective mass is above the bare mass for At — Ax^. However, as At is reduced, the effective 
mass decreases. For a 32 x 128 lattice, the effective mass is found to be around 1.07. The reason lies in the fact that 
we used the simple central difference Laplacian operator in the nonlinear part of action Si ! We expect that the perfect 
linear operator operating on a function f{x), which does not depend on t, should yield (—9^ + m)'^ f{x). However, a 
summation of the PO along the t direction does not yield the one-dimensional NA form {—d^ + m)'^, but rather has 
coefficients roughly twice that of the NA form. Therefore, it is inconsistent to simply use the central difference form 
for operator {—d^ + to). A test simulation using \/2{—d'^ + to,)na gives the value 1.36 for the effective mass, closer 
to our expectation. However, it is not clear how to interpret this and it points to the need to derive the perfect form 
for the whole action including the nonlinear part. 

For the static structure factor S{k,t = 0), shown in Figure |2|, the PO result is not very close to the bench 
mark result of NA with a 64 x 256 lattice. For large values of k, there is a contribution from the stochastic CG 
transformation. For small k values, its deviation is a result of the inaccuracy in the effective mass, which is related 

— 1/2 

to the correlation length ^ (and hence the shape of S{k)) by the relation ^ ~ "^off ■ 

It is interesting to notice that the structure factor curves obtained using different schemes and lattice sizes all cross 
at the same point around k w 1.3. 



2. The m < Regime 



In this case, there is a non-trivial fixed point in the action which corresponds to a ground state with order parameter 
values (p = zL y^m/g. Domains of opposite order parameter values compete and the dynamics is quite different from 
that with TO above 0. In our simulation, we used the same parameters as in the previous section except to, = — 1. We 



treat md 



as one term and use the massless perfect linear operator. This leads to a reduced standard error. The 



data are plotted in Figures |3| and |lj. 
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FIG. 13. Decay rates of wave modes for the Ginzburg-Landau equation. 



-1, L = 16, r = 8. 



The general shape and values of the dispersion relation are similar to those of the jti > regime. However, there is 
a marked difference between these two regimes for wave modes close to A: = 0. Here, instead of approaching a finite 
effective mass, the decay rate approaches zero, reflecting the existence of a ground state with a non-zero amplitude. 
Also due to the 'vanishing' effective mass, the shape of the structure factor is more peaked at the origin than in the 
m> regime. For modes with small k (first few nodes), S{k,uj) values have a large standard deviation. For example, 
it is about 25% for the k = An/L mode and about 9% for the k — Stt/L for NA on a 32 x 64 lattice. 
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FIG. 14. The static structure factor S{k, t = 0) for the Ginzburg-Landau equation. 



-1, L = 16, r = 8. 



When grid sizes are reduced, the dispersion relation changes shape for small k modes. The difference is significant 
with respect to the standard error. This has also been checked with increased statistics. This may be due to the 
existence of the non-trivial ground state. For m < 0, there is another length scale in the problem, namely, the interface 
width between domains with opposite signs of the ground state order parameter value. If the grid size Ax is not 
small enough, the position, and hence the dynamics, of the domain interface will not be resolved. This seems to be 
the reason why the shape of the dispersion relation for small k values changes as Ax is reduced, and it places an 
inherent physical constraint on the level of discretization one can reach. Only when this extra complication is taken 
into account can we obtain a perfect operator for this problem. Nevertheless, as shown in the figure, the perfect linear 
operator gives superior results to the NA operator for the same lattice size and computational effort (as discussed in 
the previous section). 

In summary, a direct application of the perfect linear operator gives us an improved dispersion relation for Model 
A dynamics, especially for those modes with a length scale comparable to the lattice grid size. However, a more 
extensive study is needed to fully assess the efficacy of the perfect operator. This requires improving the perfect 
operator such that it yields the correct effective mass in the m > regime and accounts for the formation of domain 
interfaces in the tti < regime. 

D. Modified 'Perfect' Operator 

As previously shown, although the perfect operator coefRcients fall off exponentially as one moves away from the 
origin, the decay rate is slow along the x direction. Therefore, an operator with a shorter range of interaction is 
desired. 

In non- linear a model [ p2| , by simply including the next-nearest-neighbors (NNN), the dispersion relation can be 
greatly improved. In that case, the NNN coefficients are obtained using a natural truncation of the perfect operator. 
Since the operator coefficients fall off quickly along both x and t axis, such a severe truncation can still lead to 
significant improvement. This is no longer true for the diffusion equation. However, we might ask, can we improve 
the NA operator by allowing for non-zero operator coefficients for more neighbors? The answer is yes. 
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TABLE V. Coefficients of the modffied perfect action operator, /i = and Aa;^/Af = 1. 
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We begin from the continuum limit constraints of equation (52). Setting /x = and keeping p{i,j) non-zero for 
(«, j) G {(0, 0), (0, 1), (1, 0), (2, 0)} (called the basic points), the conventional operator is obtained as the only solution 
to these equations. When more neighbors are included, the constraints are enforced by solving for p of the basic 
points as function of the other coefficient values. 

Using these non-basic-points coefhcients as fitting parameters, we can obtain an operator with a near perfect 
dispersion relation. If two parameters {H{1, 1) and H{2, 1)) are used to obtain a 3 x 2 operator, the average error for 
the dispersion relation is about 6%. By fitting four parameters (by also including H{3, 0) and H{3, 1)), we can obtain 
a 4 X 2 operator — called the modified perfect operator (MPO) — which yields a dispersion relation with an average 
error of 1.7% with respect to the exact result as shown in Figure Qq. The operator coefficients for m — are given in 
Table 0. 
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FIG. 15. Decay rates of wave modes for the diffusion equation using the modified perfect operator, /i = 0.25, At = Ax^. 

For the MPO, the scaling regime starts from the first time node of the two point function (i.e. S{k,t = 0)) due to 
the nearest neighbor interaction along the time direction. So long as the first two time nodes have reliable values, one 
can estimate the decay rate. This greatly loosens the precision constraint placed by the perfect operator used before. 
When the field has mass, direct fitting under modified constraints that take into account the mass cause little change 
in the coefficients. 
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FIG. 16. Decay rate of wave modes for the Ginzburg-Landau equation, m = 1, L = 16, T = 8. The modified perfect operator 
gives results comparable to that of the perfect linear operator. 



We tested the MPO in simulations of Model A dynamics. The results are comparable to that of the perfect operator 
(see Figure 16). It actually gives more accurate decay rates for wave modes at the edge of the first Brillouin zone, 
since it allows the use of the i = and t = Ai nodes to compute the decay rate, while doing this for the PO is an 
approximation. The computational effort for the MPO is drastically reduced due to the relatively short interaction 
range. 

The perfect linear operator operates on the coarse grained variable. For the modified perfect operator, the physical 
meaning of the variable it operates on is not apparent. As discussed in section I, there is a correspondence between an 
operator and a specific coarse graining scheme. For the local averaging CG scheme, or hard CG, the resulting perfect 
operator has a long interaction range. However, the range of interaction is reduced after we modify the CG scheme to 
be a soft CG scheme dependent on the parameter k. Therefore, it is reasonable to think that there is a variant of the 
standard local averaging coarse graining scheme that gives the fast decaying operator we have computed above. The 
further investigation of this point is of general interest as regards the development of an efficient numerical algorithm. 



VI. CONCLUSIONS 

The work presented in this paper is a first step towards reaping the full benefit of using renormalization group 
in the study of dynamics of spatially extended systems. We have constructed perfect representations of stochastic 
PDEs that not only integrate out the small scale degrees of freedom (in space and time) , but also develop non-local 
representations of the underlying equations that are free of lattice artifacts. We demonstrated this by computing 
the dispersion relation for elementary excitations, and comparing the results at large wavenumbers with theoretical 
expressions valid in the continuum limit. We exhibited computations for diffusion equations, and a nonlinear equation 
derived from model A dynamics, and explored different ways to truncate the non-local space-time operators generated 
by the RG. 
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In one dimension, the computational complexity was reduced by a factor of about 40 from conventional simulations, 
for the simple diffusion problem. For the nonlinear model A equation, the results were less impressive, in terms of 
computer time, because a systematic approximation scheme for the perfect action has yet to be developed. Never- 
theless, proceeding heuristically, we were still able to obtain improved results for the static structure factor and the 
decay rate of modes. Lastly, we proposed a heuristic discretization algorithm that incorporates the ideas of perfect 
operators, but also gives operators that are more local than perfect operators. 

Finding the perfect operator when nonlinear interactions are present is a non-trivial task. The form of the continuous 
action is not closed under the CG transformation and new complicated interaction terms are generated. This is a 
general property of the RG [^ . Usually progress is only possible if the problem under consideration involves a small 
parameter which can be used to keep track of the new interactions which are generated. More generally, the small 
parameter allows a systematic approximation scheme to be developed, in which there is a clear prescription as to 
which terms have to be included at a given order. If such a parameter is not available, it is usual to fall back on 
to some type of variational scheme, typically including some kind of self-consistent calculation which corresponds to 
summing sets of diagrams. Neither of these approaches have been attempted in this paper. However, we feel that the 
results which we have obtained are sufficiently encouraging that some type of systematic calculation of the perfect 
operator in nonlinear theories would turn the ideas presented in this paper into a powerful computational tool. 
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APPENDIX A: 

In this appendix, we prove various relations that are important in deriving the iterative relations for perfect 
operators. 

1. Here we list some properties of the projection matrices. 

- ~ --1 -T ~-i ~T 

We introduced 2N x N matrices R, R and their left inverses R = ^R , R — ^R , 

Rm.n = S„i^2n + <5m,2n-l, Rm,n = S,n,2n + <5m,2n-l m £ [1, 2A^], n G [1, iV]. (Al) 

Here superscript T indicates transposition. The projection matrices satisfy relations, 

--1- --1- ---1 ~--i 

R R = R R = and RR + RR = 1 . (A2) 

We define the subscripted versions of O by 

R OR = OaR OR = ObR OR = OcR OR = Od, (A3) 

where O and its subscripted versions are linear operators on the original grid and on the coarse grained grid 
respectively. 

One can prove the following formulae, 

g^Of^{ ^^ ^~^~y^^ ^"+ (^"!^^) ^~ ^OAf + Oc f (^4^ 

g^{R ^dR)f+{R 'dR)f = dDf + dBf 

fOf = 2 {F Oa I + 2/^ Oc f + F Ob /) especially f = 2 {p + p) , (A5) 
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where we assumed that the matrix O is symmetric (physicahy, this means O possess inversion symmetry), and 
therefore Od = Oq. Furthermore, if O is translational invariant with Om,n = Om+i,n+i, Oa and Ob are 
symmetric while Oc and Od are antisymmetric. This can be seen by looking at their elements, 



{OA)m,n — 02m,2n 
{OB)m.,n = 02m,2n 
idc)r. 



2{02m,2n+l + 02m,2n-l) 

5(02™, 2n+l + 02ni,2n-l) m, n = 1, ■ ■ ■ , N. 

= 2(^2m,2n+l ^ 02m,2n~l) ■ 



(A6) 



In Fourier space, (pirn) — J24'{k) exp{ikm), m E [1,2A^] and (l>{n) = ^ 0(k) exp(iKn), n E [1,^]- We have 
R = 2C.C. (R ), i? — 2C.C. (R ), where c.c is complex conjugate, and 

[ ^fc,K = \/2e'4(-zsmf4^.| +cosf4,|±^) 



The sign in -^ ± tt should be chosen so that its value lies within the interval (— tt, tt). Physically, equation (A7) 
represents a two step process: folding of the Brillouin zone by half, such that two wave modes k and fc ± tt are 
mixed, followed by a stretching back to {—it, tt). This is the corresponding process in Fourier space of the real 
space coarse graining transformation. 

Given the operator O — J20{k)\k){k\, i.e. plane wave functions form its eigenspace, the coarse grained plane 
waves |k) are also eigenvectors of Oa, Ob and Oc, 



(Oa).,k' = (Jk.k' [cos^ f 0(1 ) + sin2 f 0(| ± tt)] 
{Ob).,.' = S^.Asin^ f 0(f ) + cos^ f 0(f ± tt)] 



K, k! E {—it, it). 



(A8) 



{Oc)., 



(5,,,.(-zcosfsinf)[0(f)-0(f ±7r)] 



2. In order to determine the iterative relations of linear operators (see 3 below), we first have to prove some 
properties of the subindexed matrices. 



(i) The first set of properties are: 



{0-^)a-0c = -{0-^)c-0b 
dD-{d-^)A = -Ob -{6-^)0 

(6-1)z5-6c-i-(o-1)b-6b. 



(A9) 



To prove the first relation, we use equation ([A^): 



{o-^)a-Oc^r q-^r-r 0R = R o-^{\-rr )0R 
= r o^^or-r o^^rr or 

= R~' R- (d-i)c Ob = -{d-')c Ob ■ 



(AlO) 



We can prove the other two relations in a similar way. 
(ii) Another very useful result is: 



{d-')A = {dA-dc{dB)-'dD)-' 



(All) 



We prove this using equations (A9): 

{6-^)a ■ {Oa - dc{6B)-'dD) = {d-')A ■ Oa + {d~^)c Od 



lf^6~^RR~^dR + if^d-^RR OR 

if O-^OR 

RR^l 



(A12) 
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3. The iterative relation for the action operator (equation (19)) follows from that of p (equation (|lj)) which reads 
{pCGyi ^ LcM-^Ps\M-TlI + r{pA- PcPb'pd)-' f ^ , (A13) 



where t = I + Oc{Ob)~^ Pb^ Pd, and where we use O to denote the full dynamic evolution operator Lu,. Since 
p is symmetric, A and B subindexed mat rices are sy mmetric while C subindexed matrix is the transpose of D 
subindexed matrix. Using equations (A9) and ( All ), the second term in the above equation yields. 



T\-1AT 



r {p-')A r' = {p-')A - Oc{Ob)-' ip-')D (p-'h {o'sr'oij 

+0c(0b)-^ [{p-')b-{pb)-'] {Olr'Ol. 



(A14) 



Therefore (if for the US scheme, ignore the factor of 2), 

(/0-' = {p-')A + dc{d-')B (pb') {olr'ol 

-dciOB)-' ip-')D (p-'h {dirndl . 



(A15) 



Therefore, using the iterative relation (O^^) ^ 

\M\-1 /-M\-l t/\M\T -1 



[O )a and equations (|A9D, we have 



{o''')-'{p''')-'{o'''y 
{6-')a {p-')A (o^ -')a + {d-')A OciOB)-' {p-')B [dlr'dl (6^ -^u 
-(o-^U OciOB)-' {p-')D (o^ -')a - {6-')a {p-'h [dlr'dl (6^ -')a 
{d-')A {p-')A (o^ -')a + {d-')c {r')B (6^ -')c 

+{0-')c {p-')D (O^ -')a + {0-')a {p-')c {0^ -')c ■ 



(A16) 



On the other hand, we have 



(O-^ p-^ O^ -^)a = R 0-^{RR +RR )p-^{RR +RR )0^-^R. 



(A17) 



Expanding the above equation and comparing with equation (A16), proves equation (M 
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